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' We derive the normal form for the delay-induced Hopf bifurcation in the first-order phase-locked 

I loop with time delay by the multiple scaling method. The resulting periodic orbit is confirmed 

by numerical simulations. Further detailed numerical investigations demonstrate exemplarily that 
this system reveals a rich dynamical behavior. With phase portraits, Fourier analysis and Lyapunov 
spectra it is possible to analyze the scaling properties of the control parameter in the period-doubling 
scenario, both qualitatively and quantitatively. Within the numerical accuracy there is evidence that 
the scaling constant of the time-delayed phase-locked loop coincides with the Feigenbaum constant 
5 ~ 4.669 in one-dimensional discrete systems. 
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I. INTRODUCTION 



Many nonlinear dynamical systems in various scientific disciplines are influenced by the finite propagation time 
of signals in feedback loops. A typical physical system is provided by a laser where the output light is reflected 
and fed back to the cavity P, Q • But time delays also occur in biology due to physiological control mechanisms 
ff^ , 0, Q or in economy where the finite velocity of information processing has to be taken into account 0, @] . Further- 
■ more, realistic models in population dynamics or in ecology include the duration for the replacement of resources 0,0. 



■ All these different systems have in common that the inherent time delay may induce dynamical instabilities. 
1^ Numerous experimental and theoretical studies have demonstrated this for the emergence of oscillatory behavior, 
quasi-periodicity, chaos or intermittency. But time-delayed feedbacks can also have the opposite effect. They have 
even been devised to stabilize previously unstable stationary states or limit cycles p, 10, llj. In particular this 
method allows to control or to prevent undesirable chaotic behavior in a time-continuous way. In comparison with 
C the time-discrete method of Ott, Grebogi, and Yorke it can be easier implemented as it relies on less information 
^ . of the dynamical system. 
• 1— I . 

I The choice of a paradigmatic model system for analyzing the fundamental properties of delay-induced instabilities 
. is determined by several practical conditions. On the one hand it should be guaranteed that all observable instabilities 
■ " " ' are purely a result of time delays. On the other hand the dynamics should be governed by a simple model equation 
and allow for a quantitative comparison between theory and experiment. These conditions are fulfilled, for instance, 
by the Mackey-Glass model Q which describes quite successfully the anomalies in the regeneration of white 
blood cells. However the underlying nonlinear scalar delay differential equation necessitates three effective control 
parameters to account for the experimental data. 

In this article we report about analytical and numerical investigations of the rich dynamical behavior of another 
model system which was proposed some time ago by Wischert et al. in Ref. I* represents the first-order 

phase-locked loop (PLL) with time delay which synchronizes the phases of two oscillators. In comparison with the 
Mackey-Glass model, this system has the advantage that it involves only a single effective control parameter instead 
of three. Additionally, it can be realized electronically under well-defined conditions. Furthermore, an extension of 
the PLL with time delay is capable of describing sensibly physiological control experiments [l^ . 
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The experimental set-up for the electronic system of a first-order phase-locked loop (PLL) is shown in Fig. 3 
of Ref. In many applications the PLL serves for synchronizing the phases of a reference oscillator and a 

voltage-controlled oscillator (VCO). Thereby the frequency of the output signal of the VCO depends linearly on the 
input signal. The output signals of both oscillators are multiplied by the aid of a mixer. The induced high-frequency 
components are then eliminated by a low-pass filter. The resulting signal is fed back to the input of the VCO. A 
delay line between the VCO and the mixer, implemented analogously or numerically, induces the time delay t > 0. 

The dynamical variable of interest is the phase difference q{t) between both incoming signals of the mixer (compare 
with Fig. 3 of Ref. 0). Under quite simple assumptions it becomes possible to derive a nonlinear scalar delay 
differential equation for this phase difference [Tsf : 

j^q{t) = -K sin[g(< - r)] . (1) 

The parameter K >Q denotes the so-called open loop gain of the PLL. Performing an appropriate scaling of the time 
converts the PLL equation to its standard form 

j^q{t) = -R sm[q{t - 1)] , (2) 

where the two parameters t, K are reduced to one effective control parameter 

R^Kt. (3) 

Thus varying the delay time t corresponds to changing the control parameter R. In this paper we analyze both 
analytically and numerically the rich structure of delay-induced instabilities of the PLL equation In Section^ we 
derive the normal form of the Hopf bifurcation by applying the multiple scaling method. The emerging periodic orbit 
is confirmed in Section lllll bv numerical simulations. In Section llVI we study in detail the period-doubling scenario 
beyond the Hopf bifurcation with phase portraits, Fourier analysis and Lyapunov spectra. 



II. MULTIPLE SCALING METHOD 



Applying the synergetic system analysis it was shown in Ref. [U] that a Hopf bifurcation occurs in the PLL equation 
(O at the critical value Rr = 7r/2. In the following we rederive the normal form of this Hopf bifurcation by using the 
multiple scaling method |25j. It represents a systematic technical procedure to deduce the normal form by using an 
ansatz how the respective quantities depend on the smallness parameter 

^^?^_Jh ^ i? = i?,(l + e). (4) 

Althou gh the multiple scaling method has been originally developed for ordinary or partial differential equations 
|26L I27L |28L I29I I30II . it can be also applied to delay differential equations (see, for instance, the treatment in Ref. 31j). 



We start with discussing some properties of the Hopf bifurcation. At first, we mention that the amplitude of the 
emerging periodic solution has a characteristic e-'^/^-dependence from the smallness parameter e, as can be deduced 
from the linear stability analysis already presented in Ref. Furthermore, the trajectory approaches the limit 

cycle slowly near the instability due to the phenomenon of critical slowing down. Thus the oscillatory solution q{t) of 
the PLL equation ||2Jl is based on two different time scales. The fast time scale is provided by the period T = 2n/Vl 
of the oscillatory solution, whereas the slow one characterizes the amplitude dynamics in the transient regime. Both 
time scales are separated by a factor of the smallness parameter e, as follows again from the linear stability analysis 
[l3l| . These considerations lead to the following ansatz for the oscillatory solution after the Hopf bifurcation: 

q(t) - <Zstat + e"^ [qUt')e^''' + 9o"(i')e-*"*] + e [q+ {t')e^^''' + q^{t') + q- {t')e-^^^'] 

+£3/2 [q+{t')e''^^ + g+(t')e*''* + q^{t')e-'''' + (t')e~''"*] + O (e^) . (5) 

Here the first time scale t and the second time scale t' are related via 



t' = et, 



(6) 
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where the smallness parameter e denotes the deviation from the bifurcation point according to Q , and the respective 
amphtudes have the properties 

9^(0 = 9,^(i')*, A: = 0,2, 3, 4; gi(i') - ■ (7) 

Now we insert our ansatz (0) in the PLL equation By doing so, we have to take into account that the slowly 
varying amplitudes qf{t') have the time derivative 

l,*(,,.4,f(f) (8) 

and that their time delay results in the expansion 

iHt' - = it it') - e ^qtit') + Oie') . (9) 

The strategy is then to compare all those terms with each other which have the same power in the smallness parameter 
e. Thereby we have to guarantee that in each order the respective Fourier coefficients compensate each other: 

• In the lowest order e° we have only the frequency which leads to the equality 

i?c Sin(feat) =0. (10) 

This fixes the stationary state to be 

gstat = Z^, Z = 0,±1,±2,... . (11) 

In the following we restrict ourselves to consider the reference state ^st^t = as the other one qH^^^ = tt turns 
out to be unstable for all values of the control parameter R. 

• The order e^/^ contains only the frequency ±fl with the condition 

±inq^{t') = -R,e^'''qt{t'). (12) 

As the amplitudes q^it') should not vanish, we conclude 

- i?ce=F''" Ti^^O- (13) 

This condition coincides with the transcendental characteristic equation (98) of the linear stability analysis of 
Ref. 01 if the eigenvalues A at the instability are identified according to A = ifl. The real part of Eq. H13|l 

cos 17 = (14) 

leads to the frequency 

f7=|, (15) 
whereas the imaginary part results in the critical value of the control parameter: 

Rc^l- (16) 

• The order e involves two frequency components. The Fourier coefficients of the frequency immediately lead to 

9i(i')-0, (17) 

and for the frequency ±217 we obtain 

± 2inqf{t') = -R,e^^'^\tit') , (18) 
which reduces due to the characteristic equation (|13|1 to 

Q2±(t')-0. (19) 
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Also the order e'^/^ consists of two frequency components. For the frequency ±51 we read off 



(1 - i?,eT'^") -^q^it') = {~R,e^^^ T ^^^) qHt') - ^Rce^''' [2qtit') - qHt'f^it') 



(20) 



The factor in front of q^{t') vanishes because of the characteristic equation Thus the functions q^{t') 

are not determined in this order, they only follow from the next order and the frequency ±fl. Taking into 
account the characteristic equation we yield from (|20l) the normal form of the order parameter equation: 



There the parameters and are defined by 

1 + iRc' 2 ■ 



A' 



(21) 



(22) 



Note that the normal form (|21|1 of the multiple scaling method and the normal form of the synergetic system 
analysis in Ref. ^3] do not coincide, however, they can be mapped into each other using an appropriate coordinate 
transformation j25(. Correspondingly, we obtain for the frequency ±357 



which reduces due to H13|l . (|15|l and to 



qfit" 



(23) 



(24) 



As the functions qf{t') and qf(t') are of the order e'^/^ according to the ansatz 0, they are irrelevant for the 
order e in which we are interested. 



It remains to solve the order parameter equation H21() . H22|l by using polar coordinates 

q^{t')^R{t')e^'^'^''K 
The resulting stationary solution turns out to be 

R{t') = V2 + 0{e^) 



together with the phase 



Here the frequency turns out to be 



ip{t) = n{e)t + ifio . 



(25) 
(26) 
(27) 
(28) 



Thus we conclude from lfTB|) , ifTTjl , ((T^ , (|5SJ , if^ , and (|2Hl that the oscillatory solution after the Hopf bifurcation 
reads 

q{t) = co(e) + ci(e) cos [ip{t) + V'l] + £2(2) cos [2^{t) + + C (e^/^) , (29) 

where the respective coefficients read 

co(e) = + O(e2) , ciie) ^ Vs^ + O (e^^^^ , €2(2) - + O (e^) . (30) 

It coincides with the result of the synergetic system analysis in Ref. up to the order e. Although we restrict 
ourselves to this order, the systematics of the multiple scaling method is obvious, thus an extension of the ansatz © 
to higher orders is straight-forward, but the calculation would become quite cumbersome. 
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TABLE L Comparing the analytical and the numerical values for the frequency 0.{e) and the Fourier coefficients co(e), ci(e), 
C2(e) of the oscillatory solution of the PLL equation after the Hopf bifurcation. 



III. NUMERICAL VERIFICATION 



In order to numerically verify our analytical result, we integrated the underlying PLL equation (0). By doing 
so, we varied the control parameter R in the vicinity of the instability = 7r/2 in such a way that the smallness 
parameter e = {R — Rc)/Rc took 200 equidistant values between 10"^ and 10"^. We used a Runge-Kutta-Verner 
method of the IMSL library as an integration routine and performed a linear interpolation between the respective 
values. In particular in the immediate vicinity of the instability the phenomenon of critical slowing down led to a 
transient behavior. To exclude this, we iterated the discretized delay differential equation for each value of the control 
parameter at least 10^ times. Afterwards we calculated the power spectrum with a complex FFT so that the basic 
frequency SI of the oscillatory solution could be determined with high resolution. Then we performed a real FFT with 
the period T = 27r/r2 of the simulated periodic signal q{t) = q{t + T): 

oo 

q{t) = ^ + ^ [flfe cos {km) + hk sin {kVLt)] . (31) 

k=l 

The Fourier coefhcients follow from integrations with respect to one period T — 27r/Sl: 

T T 

ak^ ^ J dtf{t) COS (knt) , fc^O,l,...,oo; bk ^ ^ J dt f {t) sin (knt) , /c = 1, . . . , cx) . (32) 



From H31|l follows then the spectral representation 

oo 

q{t) = Co +'^Ck COS {kftt + (f>k) (33) 



fc=i 



with the quantities 



Co = ^ , Ck = Jal + bl, (^fc = -arctan — , fc = 1, . . . , oo . (34) 



Thus our analytical result (|27|l . (|29() can be interpreted as the first terms within a spectral representation 133|l . 
where the frequency H. = 27r/T and the Fourier coefficients Cq, Ci, C2 are given by (|28|l and (|30|l . Analyzing the 
Hopf bifurcation with a FFT, we numerically determinded fi, Cq, Ci, C2 as a function of the smallness parameter e. 
Comparing the respective numerical and analytical results, we observe some deviations for small and for large values 
of the smallness parameters e. The former are due to the phenomenon of critical slowing down, i.e. the system stays 
longer in the transient state when the instability is approached, and the latter arise from the neglected higher order 
corrections in the analytical approach. Therefore we restricted our numerical analysis to the intermediate interval 
[10~^, 10^^] of the smallness parameter e. 
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In Tab. ^ we see that the analytical and numerical determined quantities agree quantitatively very well. Thus our 
weakly non-linear analysis for the delay-induced Hopf bifurcation in the PLL equation is numerically verified up to 
e ~ 10~^. However, this successfully tests only the order parameter concept for delay systems, as the lowest nonlinear 
term in the scalar delay differential equation of the PLL (j^J is a cubic one. Therefore we analyzed the Wright equation 
[32| with a quadratic nonlinearity in a separate publication , where we could successfully test not only the order 
parameter concept but also the slaving principle, i.e. the influence of the center manifold on the order parameter 
equations. Thus we demonstrated with the Wright equation the validity of the circular causality chain of synergetics 
for the Hopf bifurcation of a delay differential equation. 



In the following we summarize various simulations which have been performed for the delay differential equation 
(0) of a PLL with some initial function q{t) for — 1 < t < [25ll3^. In order to check the quality of the numerical 
results, standard integration routines of the Runge-Kutta type have been applied with different discretizations by 
adequately taking into account the delay effects. 

Firstly, it turns out that the trajectory q{t) is restricted for all times to the interval [—tt, +7r] if the effective 
control parameter R is increased from to about 4.9. If the transient behavior has decayed, the resulting asymptotic 
dynamics could be classified as follows. For < R < n /2 there exist two stationary states, a stable one gi = and 
an unstable one (72 = tt. At i? = 7r/2 a super-stable Hopf bifurcation occurs where the previously stable stationary 
state (7sta^t = becomes unstable and a new stable limit cycle emerges In the range tt/2 < R ^ 3.77 this 

oscillatory solution shows a conspicuous point symmetry with respect to the origin of the phase portrait in Fig. 
This symmetry is broken at i? « 3.77 as the limit cycle splits into two coexisting limit cycles [s^- They are depicted 
in Fig.^ as the asymptotic dynamics of the initial functions q{t) = ±1 for —I < t < 0, respectively. Both coexisting 
limit cycles are symmetric to each other concerning the point symmetry with respect to the origin and remain stable 
up to i? w 4.9. Note that the instability at i? w 3.77 was not detected during the initial investigations in Ref. |0 as 
there only Fourier spectra were analyzed. 

Increasing the effective control parameter leads to a further bifurcation at i? w 4.105. Fig. ^ illustrates that a 
new limit cycle emerges with the initial function q{t) = 2 for — 1 < t < which coexists for 4.105 ^ R ^ 4.11 with 
the two other oscillating solutions. Also this new limit cycle exhibits a point symmetry with respect to the origin 
of the phase portrait. At i? « 4.11 this limit cycle splits into two new oscillating solutions with the initial functions 
q{t) = ±2 for — 1 < i < 0, so that the point symmetry is again destroyed (compare Fig. ^). It turns out that 
both of them pass through a separate period-doubling scenario for 4.11 ^ R ^ 4.175. This is shown qualitatively by 
the Power spectra (see Fig. |2l for the first three period-doublings. Each of these bifurcations leads to a subsequent 
sub-harmonic and to corresponding higher combination frequencies. The bifurcation diagram in Fig.|21is an overview 
over this period-doubling scenario. It was obtained by Poincare sections of the trajectories using the software package 
AnT 4.669 l4ll,l42], whereby the Poincare conditions were q{t) = and q{t) e [1,2]. 

In order to analyze a period-doubling scenario more quantitatively, it is advantageous to determine the Lyapunov 
exponents of the underlying dynamics. In our case this necessitates to use the common concept of Lyapunov 
exponents [36| and to extend it to delay differential equations 25. .37j . Thereby we have to take into account 
that their numerical integration is based on a discretization procedure. As a consequence, the determination of 
Lyapunov exponents for time-delayed dynamical systems is reducible to the calculation of Lyapunov exponents for a 
high-dimensional time-discrete mapping. 

Fig.^ shows the two largest Lyapunov exponents within and above the period-doubling scenario of the PLL. The 
enlargement of Fig. 0Jd clearly reveals the self-similarity of the spikes and the characteristic scaling properties. One 
of both Lyapunov exponents always vanishes due to the moving reference frame. The zeros of the second Lyapunov 
exponent coincide with the critical values of the effective control parameters i?o, i?2, . . • where a period-doubling 
occurs. Tabled lists the first bifurcation points and the scaling constants 



of the effective control parameter. Within the numerical accuracy there is evidence that the scaling constants Sn 
converge to the Feigenbaum constant S ~ 4.669 as in the case of one-dimensional time-discrete systems [38L Is^. If 



IV. FURTHER NUMERICAL RESULTS 



Sn = 



R. 




(35) 



R. 
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TABLE II: Bifurcation points of the period-doubling scenario. 



we assume this to be true then we can estimate the end of the period-doubhng scenario 

= ^" (36) 

from the bifurcation points in Tab. Hll The finding i?oo ~ 4.173961 agrees quite well with the enlarged Lyapunov 
spectrum in Fig. ^p. 
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FIG. 1: Phase portraits depicting several limit cycles: a) R — 3.5, b) 7? = 4.1, c) R — 4.105 d) R — 4.11. 



As typical for the period doubling route to chaos, there exist not only the period doubling scenario below the 
critical value, that is for R < Roc, which ends at the critical value i?oo with the emergence of a Feigenbaum attractor, 
but also the band-merging scenario with periodic windows above the critical value, that is for R > i?oo- As an 
example the periodic behavior in the window 4.2095 ^ R ^ 4.215 was analyzed more carefully. The phase portrait of 
Fig. inii and the Power spectrum of Fig.jSJj show that at i? « 4.2095 a limit cycle of period 3 emerges. At i? « 4.2115 
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FIG. 2: Power spectra indicating period-doublings: a) R = 4.157, b) i? = 4.165, c) R = 4.1725, d) R ^ 4.17375. 




FIG. 3: Bifurcation diagram obtained by a Poincare section which is defined by the conditions q{t) = and q{t) £ [1, 2]. 



it starts to pass through a period-doubling scenario (see Fig. [S];). Finally, at i? « 4.213 a chaotic attractor emerges 
whose power spectrum in Fig. (Sji clearly reveals the structure of the limit cycle of period 3. This chaotic regime 
ends at i? « 4.2405 when a global bifurcation or a transition to transient chaos occurs. For 4.2405 ^ R ^ 4.85 phase 
portraits and Power spectra show that only those limit cycles coexist which emerged at i? w 3.77. At i? w 4.85 two 
new limit cycles of period 2 are generated which coexist for 4.85 ^ R ^ 4.90. 

For < R ^ 4.90 the dynamics has the characteristic property that the state space F is divided in separate 
intervals [(m — l)7r, (m + l)7r] with m = 0, 1, 2, . . .. In each of these intervals occurs the dynamical scenario which 
has been described so far. At i? « 4.90 it happens for the first time that previously separated intervals are linked 
together, so that a new dynamical behavior becomes possible. For 4.90 ^ R ^ 5.30 the Figs. and show 
that there exist, for instance, limit cycles of period 2 in different intervals although the constant initial function 
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4.15 4.175 4.2 4.225 4.25 4.1737 4.1738 4.1739 4.174 

7? R 



FIG. 4: Lyapunov spectra illustrating self-similarity. Shown are the two largest Lyapunov exponents Ai and A2. 

q{t) — —2 (dashed-dotted) —1 (dashed), 1 (dotted), and 2 (solid) was chosen in the interval [— tt, +7r]. Thus the 
transient dynamics occurs in different intervals, whereas the asymptotic dynamics is restricted to one of these intervals. 

At i? ~ 5.30 occurs another instability to chaotic behavior. However, now the chaotic dynamics is no longer 
restricted to one of the above mentioned intervals, but it relates the previously separated intervals. For a certain time 
the system dynamics remains restricted to one of these intervals and moves then to the next interval (see Figs.lHt and 
EJi) . Thereby the time duration of the system within one interval differs from interval to interval. Such a dynamics 
is called phase slipping or cycle slipping All numerical investigations above R « 5.30 show that only the phase 
slipping behavior is stable. Analyzing the Lyapunov spectrum Ai > A2 > • • ■, the Lyapunov dimension 

j 

is found to increase linearly with the control parameter R (compare Fig. [Tj). Note that also other time-delayed 
dynamical systems possess chaotic attractors where the envelope of the Lyapunov dimension Dl is proportional to 
the time delay r |25l ISTl l40l | . To our knowledge a theoretical explanation for this universal phenomenon, which could 
predict the system specific slopes from the respective delay differential equations, is still lacking. 

V. CONCLUSION 

Here we have demonstrated by the example of the phase-locked loop with time delay that an adequate combination 
of different analytical and numerical investigation methods reveals different aspects of the rich dynamical behavior 
in time-delayed nonlinear systems. The multiple scaling method allows to derive the normal form for the Hopf 
bifurcation. Phase portraits resulting from different initial functions are capable of detecting the splitting of a limit 
cycle indicating thereby symmetry-breaking bifurcations. A period-doubling is qualitatively indicated in the power 
spectrum, whereas the Lyapunov spectrum allows more quantitative statements. In this way we found within the 
numerical accuracy evidence that the period-doubling scenario in the phase-locked loop with time delay is governed 
by the Feigenbaum constant S « 4.669. 
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